Phase transition in the three dimensional Heisenberg spin glass: 

Finite-size scaling analysis 
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We have investigated the phase transition in the Heisenberg spin glass using massive numerical 
simulations to study larger sizes, 48 , than have been attempted before at a spin glass phase tran- 
sition. A finite-size scaling analysis indicates that the data is compatible with the most economical 
scenario: a common transition temperature for spins and chiralities. 
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I. INTRODUCTION 

As a result of extensive and careful numerical 
studies^^i 3 - there is now compelling evidence for a fi- 
nite temperature phase transition in the Ising spin glass 
in three dimensions. However, the situation for the 
Heisenberg spin glass, in which the spins are classical 3- 
component vectors, is still controversial. The Heisenberg 
spin glass is a suitable first model to describe experi- 
mental systems with weak anisotropy, such as dilute Mn 
atoms in Cu which is a well studied spin glass system, see 
e.g. Ref. H Kawamura^ proposed that the spin glass 
transition only occurs at Tgc = and that a chiral glass 
transition occurs at a finite temperature Tcg ■ Chiralities 
are Ising-like variables which describe the handedness of 
the non-collinear spin structure. This scenario requires 
that spins and chiralities decouple at long length scales. 
However simulations 7 -^^ subsequently found evidence for 
a finite Tsg, though corrections to the leading finite-size 
scaling behavior seem larger than in the Ising case4 Re- 
cently, Viet and Kawamur a 10 ' 11 who did a similar anal- 
ysis to that of Refs. 3^ and used the same range of sizes 
(L < 32, where L is the linear size of, the system), con- 
cluded that Tgc is indeed finite, but is less than Tqg 
which still implies spin-chirality decoupling. 

In view of this controversy over the nature of the tran- 
sition in the three-dimensional Heisenberg spin glass, 
which is of great importance for the understanding of spin 
glasses, we have undertaken a massive set of simulations 
to study even larger sizes^ N = L 3 where L < 48. Our 
conclusion is that the data is consistent with a common 
transition temperature for spins and chiralities, though, 
of course, numerics can never prove that they are exactly 
equal. 

The paper is organized as follows. In Sect. [II] we de- 
fine the model and the observables. Finite size scaling, 
which is central in our analysis, is recalled in Sect. IIIII 
Simulation details are in Scct. lIVl while our equilibration 
tests are addressed in Sect. IVl We find that a uniform 



allocation of computational resources is inefficient (equi- 
libration is much harder to achieve for some particular 
samples). The numerical results are in Sect. IVI1 while 
our conclusions are presented in Sect. IVTT1 



II. MODEL AND OBSERVABLES 

We use the standard Edwards-Anderson spin glass 
model on a cubic lattice 



T~i — S ' Jij Si ■ S 

(i,j) 



J ■ 



(1) 



where the Si are 3-component classical vectors of unit 
length at the sites of a simple cubic lattice, and the Jij 
are nearest neighbor interactions with a Gaussian dis- 
tribution with zero mean and standard deviation unity. 
Periodic boundary conditions are applied. 

The spin glass order parameter is q? v = S" , 
where "(1)" and "(2)" are two identical copies of the 
system (same interactions), and v and // are spin compo- 
nents. Its Fourier transform at wave vector fe, is denoted 
by q^(k). 

For the Heisenberg spin-glass, Kawamura 5 - defines the 



chirality from spins on a line: 



i+fi 



where here fx refers to a direction on the lattice. The 
related chiral spin-glass order parameter is q GG i = 



its Fourier transform being q GG (k). 



The wave vector dependent susceptibilities are com- 
puted from the two order parameters: 

XsgW = iV]T[(|<r(fc)| 2 >] a v, (2) 

XcG (k) = N[(\q£ G (k)f)] av , (3) 

where (• ■ • ) denotes a thermal average and [• • ■ ] av denotes 
an average over disorder. 
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The susceptibilities yield the second-moment finite- 
lattice estimator of the correlation lengthfi i 13 i 14 i 15 



i 



/ x(Q) 

2sin(7r/i) V x( k min) 



(4) 



with fc m in = (2ir/L,0, 0). The spin and chiraU^ correla- 
tion lengths are denoted by £sg,l and £cg,l respectively. 

We also consider the spin and chiral Binder ratios de- 
fined by 
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where g 2 = £^ W(k = Q)]\ q CG = ^ [^ G (fc = 0)] 5 



III. FINITE-SIZE SCALING 

Finite size scaling is a most useful data analysis 
method, that exploits finite size effects where they are 
largest (at criticaiity) to gather information on the infi- 
nite system, see e.g. Ref. [la . 

Finite size scaling takes the form of an asymptotic ex- 
pansion on the system size, L. To leading order, for 
a quantity O diverging in the thermodynamic limit as 
O oc IT - T c \ xo , it takes the form 



0(L,T) = L xa ' v f (l}' v {T -T c )) , 



(6) 



where / is an analytic function of its argument. In par- 
ticular, since ^sg.l/L is dimensionless it has the finite 
size scaling forma l 17 ' 18 



£sg,l 



X (l 1/v {T -T E 



SG, 



(7) 



where v is the correlation length exponent. There are 
similar expressions for £cg,l/T, and also for the Binder 
ratios in Eq. ((5]) since these too are dimensionless. 

We shall see that corrections to scaling are quite large, 
even for the large sizes that we study, and so we need 
to consider corrections to the asymptotic scaling form 
in Eq. ([7]). To investigate this we determine the in- 
tersection temperatures Tg G (L, sL), where the data for 
£sg,l /L agree for sizes L and sL, i.e. 



$SG,L _ jsG.sL 

L ~ sL 



(8) 



with an analogous expression for the chiral data. 

Whereas Eq. ([7]) predicts that all the T£ G (L,sL) are 
equal to the spin glass transition temperature Tgc, 
when one includes the leading corrections to scaling the 
Tg G (L, sL) are given byi^i 8 - 



T* G (L,sL)-T SG 



A sg L 



(9) 



Here, uj is the exponent for the leading correction to scal- 
ing while the amplitude is 



A SG = A SG 



1 - S~ u 
S 1 /" - 1 ' 



(10) 



with A$q a (non-universal) constant. In practice, we 
do not have enough information to determine the s de- 
pendence in Eq. ([9]), so we take the Ag G to be separate 
constants for each value of s that we use (s = 2 and 3/2). 

In fact we may combine Eqs. © and ©, to ob- 
tain a modern form of Nightingale's phenomenological 
renormalization^ the so called quotient method^ 



Q(sL,T^ G (L, S L)) 
0(L,T* G (L,sL)) 



„Xo/vsG 



l+An,L- 



(11) 



where the dots stand for higher order scaling corrections. 
Were Tqg and Tsg to be different, a similar expression 
would hold for quantities diverging at T GG . In particu- 
lar, one may use Eq. (| 1 1|) with temperature derivatives 
(to obtain 1 + l/i/) or with the susceptibilities at zero 
wave number (to obtain 2 — ?y, where rj is the anomalous 
dimension) . 



IV. SIMULATION DETAILS 

Simulations are run in parallel on Nt processors at 
Nt different temperatures using the parallel tempering 2 - 
(PT) method to speed up equilibration at low-T, see Ta- 
ble |T] for the parameters. 

As discussed in other work ; 8 ' 9 ' 10 three types of moves 
are performed: (i) "overrelaxation" (OR) sweeps (which 
do not change the energy), (ii) "heat-bath" (HB) sweeps 
(which do change the energy), and (iii) parallel tempering 
(PT) sweeps in which the spin configurations at neighbor- 
ing temperatures are swapped with a probability which 
satisfies the detailed balance condition. It is important to 
include OR sweeps because not only is the code for them 
much simpler (and hence faster) than that for the HB 
sweeps, but also because OR moves are very efficient&2i2i 
in equilibrating the system, so many fewer sweeps are re- 
quired than in a simulation with only HB and PT moves. 
Nonetheless, a fraction of the moves must be HB in or- 
der to change the energy. Because of the PT moves, the 
temperature of a set of spins (a "copy" ) is not fixed but 
does a random walk between the minimum and maximum 
temperature in the set. In this way, each copy of spins 
can visit different "valleys" in configuration space with 
the correct statistical weight, even at low temperature. 

The Nt temperatures were arranged in a geometric 
progression between T m ; n = 0.12 and T max = 0.19. We 
do 1 HB sweep followed by 5L/4 OR sweeps and then 100 
PT sweeps. We found a net CPU gain by doing a num- 
ber of OR sweeps between HB sweeps which is somewhat 
greater than L, perhaps because this transfers a fluctua- 
tion right across the system. We do a large number (100) 
of PT sweeps following right after each other, because the 
PT sweeps are very inexpensive in CPU time. 
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TABLE I: Parameters of the simulations. At is the number 
of temperatures. For the larger sizes, the number of sweeps 
varied from sample to sample. We show values for N£££ ep and 
Af™ p , the minimum and maximum number of overrelaxation 
sweeps. 



L 


/V mm 

- 1 ' sweep 


Arm ax 
'sweep 


N T 


N 

1 ' samp 


8 


5.0 x 10 b 


5.0 x 10 b 


5 


984 


12 


7.5 x 10 5 


7.5 x 10 5 


9 


984 


16 


1.0 x 10 6 


1.0 x 10 6 


15 


984 


24 


1.5 x 10 6 


1.2 x 10 7 


27 


984 


32 


4.0 x 10 6 


1.2 x 10 8 


43 


984 


48 


6.0 x 10 7 


6.0 x 10 s 


79 


164 




t io' 



0.01 







_ N _ ^sample 
sweep 


A n = A min 

sweep 





0.1 



t/N 



FIG. 1: (color online) Top: Parallel tempering autocorrela- 
tion function,— as computed for three representative L — 48 
samples. Here easy means that after N^^ ep MC steps, see 
Table [TJ the equilibration criterion was met (42% of sam- 
ples), while medium samples (34%) required up to 2N££" ep 
MC steps. Bottom: The quantity A in Eq. (|12|l as a func- 
tion of MC time, both in units of N£^£ ep (red triangles) and 
in units of the maximum number of sweeps for each sample 
(blue circles). For the latter, note that the data is computed 
at different times for different samples. 



V. EQUILIBRATION 



which is valid for a Gaussian bond distribution. Here 

U = S.-S, . (13) 

g, = (l/Ay^T^.S,) 2 , (14) 

q s = (l/iV b )]r<(S 4 .S,) 2 ), (15) 



in which U is the thermally averaged energy per spin, q\ 
is called the "link overlap", A^ = (z/2)N is the num- 
ber of nearest neighbor bonds, and z (= 6 here) is the 
lattice coordination number. Both U and q s , being a sin- 
gle thermal average, come close to equilibrium relatively 
quickly as the number of MC sweeps increases. However, 
qi involves a double thermal average, which is determined 
from two separate copies initialized with random spin 
configurations, and hence is initially very small. As the 
simulation proceeds, q\ increases towards its equilibrium 
value, so A in Eq. (p~2|) will initially be positive but will 
become zero (and stay zero) when equilibrium is reached. 
Data is shown in Fig. Q] for L = 48 (the largest size) at 
T = 0.12 (the lowest temperature). 

Eq. (|12p also provides a control variat^^ to reduce sta- 
tistical errors in £sg.l and £cg.l- The key is in the strong 
statistical correlations between the Monte Carlo estima- 
tor for A and those for the susceptibilities. Since we 
happen to know that A = 0, reduced- variance estimators 
for the susceptibilities are obtained straightforwardly (see 
Ref. for details). In practice, this method halves the 
errors for £cg,sg at T = 0.12 (however, for T > 0.14 the 
gain is less than 10%). 

Some samples are harder to equilibrate than others so, 
ideally, we should spend more MC sweeps on the "hard" 
ones than on "easy" ones. The key to classifying sam- 
ples in a PT simulation is to consider the dynamics of 
the temperature random-walk. In "hard" samples the T 
random-walk is slower (a copy trapped in a deep valley 
needs a longer time to wander to a T high enough to es- 
cape). We use correlation functions and autocorrelation 
times to formalize this idea, see Fig. [T] and the comments 
in Ref. |H 

For each sample, we impose a minimum number of 
sweeps, TablelH then keep simulating until the total num- 
ber of MC iterations exceed 9 autocorrelation times. For 
L = 48, the average number of MC iterations per sam- 
ple was 1.8 times the minimum. Figure [T] shows that 
the data equilibrates more convincingly by running the 
"hard" samples for longer than the "easy" samples. 



We do several tests to ensure equilibration. Firstly, we 
require that data satisfy the relation^i 2 ^ 



A = 



T 



2 -u 

z 



= 0, 



(12) 



VI. RESULTS 

We now present our results. Figure [2] shows data for 
the spin glass and chiral glass correlation lengths divided 
by L. The resulting intersection temperatures, obtained 
from a jackknife analysis, are shown in Table ITT1 
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FIG. 2: (color online) Data for the spin glass and chiral glass 
correlation lengths divided by system size. For L — > oo the 
data should intersect at the transition temperature. Here, 
the data does not show a common intersection temperature, 
indicating that there are strong corrections to scaling for the 
range of sizes studied. 



TABLE II: Table of intersection temperatures T*(L,sL) for 
the spin glass and chiral glass correlation length data pre- 
sented in Fig. [2] Also shown are estimates for the exponents v 
and 2 — r\ for the case of s = 2, using the quotient method , 15 ' 18 
see Eq. (fTTj) . The operators used are <9t£sg,l , <9t£cg,l xsg 
and xcg which have scaling exponents 1 + 1/^sg, 1 4- 1/^cg, 
2 — ?7sg and 2 — i)cg respectively. Spin (chiral) exponents are 
computed from data at T$ G (L, sL) (Tq G (L,sL)). 

L sL Tg G Tg G z^sg ^cg 2 - wsg 2 - 77CG 

8 16 0.158(1) 0.156(1) 1.01(2) 1.34(5) 1.99(1) 0.72(2) 

12 24 0.142(2) 0.150(1) 1.35(5) 1.51(6) 2.08(1) 0.96(3) 

16 32 0.136(1) 0.147(1) 1.50(7) 1.46(6) 2.14(1) 1.11(3) 

24 48 0.133(2) 0.142(1) 1.49(13) 1.30(8) 2.19(2) 1.44(4) 

8 12 0.164(2) 0.157(2) 

16 24 0.135(2) 0.147(2) 

32 48 0.130(3) 0.138(2) 



Since Eq. ((9]) holds only asymptotically, for large L, 
it is necessary to decide on the smallest size T m i n to be 
included in the analysis. We consider first the five pairs 
of sizes with L nlin = 12, see Table [Til Fitting spin and 

chiral data separately there are 4 parameters for each: 

(2) 

Tcg.sGj the exponent uj + and amplitudes A C q sg 

and A y C Q g G for the s = 2 and s = 3/2 size ratios. We de- 
termine the best fit parameters, and estimate the quality 
of the fit frorn^ the value of % 2 . Fitting the spin data 
to Eq. ©, gives T SG = 0.1291°°°]?, which is compatible 
with Viet and Kawamura's result of 0.120(6). However, 
in the chiral sector x 2 as a function of Tcg does not have 



a local minimum with Tcg > 0, so sublcading scaling 
corrections are sizable for chiralities and T m i n = 12. 

Hence we have also performed an analysis with a larger 
value, L m i n = 16. Unfortunately, we only have data for 
four pairs of sizes, and still four parameters to be fitted 
if we fit the spin and chiral data separately. Since the 
number of points is equal to the number of parameters 
we do not gain useful information. However, if we as- 
sume a common transition temperature and do a joint fit 
we have 8 data points, and 6 parameters (1 transition 
temperature, one exponent, and 4 amplitudes). The re- 
sulting fit gives T c (= T SG = T CG ) = 0.120±°' ° 10 „ with a 
X 2 per degree of freedom of 0.029 so the fit is good. The 
error bar on T c is very large on the low-T side but if wc 
assume that the exponent \jv + lu in Eq. © is greater 
than 0.5, plausible given the values for v in Table HT1 wc 
find T c = O.^O^qq", much more tightly constrained. 

In Fig. [3] we show data for the spin glass and chiral 
glass Binder ratios defined in Eq. ([5]). Our definition of 
gee differs from that of Kawamura and Vietji - and our 
results have an even more pronounced negative dip. In- 
terestingly, we find that the results for gsG olso become 
negative at the largest sizes. Hence, the apparent vanish- 
ing of gsG near Tcg , a strong argument for spin-chirality 
decoupling^ is an artifact caused by the lattice sizes be- 
ing too small and the temperatures too high. Our inter- 
pretation of the Binder parameter data is that there is 
negative dip in both channels, and the minimum of this 
dip approaches the transition temperature as L grows. 
The chiral-dip approaches Tcg from high temperatures, 
while the spin-dip approaches Tsg from low tempera- 
tures, where plausibly Tsg = ?CG- However, much larger 
sizes would be needed to confirm this hypothesis. 

By studying sizes L < 32 Viet and Kawamur a 10 ' 11 
find T SG = 0.120(6) and T CG = 0.145(4). Since the 
difference is about 3.5 times the errors they argue that 
Tcg > Tqq. However, their value for Tcg is actually 
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higher than our intersection temperatures for L = 48 
shown in Table [IIJ an d so seems to us to be too high. 
Also they estimate the transition temperatures from 
T*(L,sL) — Tqg = const. /L av , where L av is the aver- 
age of L and sL, rather than Eq. ([9]). In other words 
they replace the exponent w + 1/y by 1, and the s de- 
pendence in Eq. ([9]) by 2/(1 + s). At the very least, we 
argue that these replacements lead to an underestimate 
of the error bars. Hence, we do not feel that the results of 
Viet and Kawamura contradict our conclusion that the 
data is consistent with the spin and chiral glass transition 
temperatures being equal. 

VII. CONCLUSIONS 

In summary, our low-temperature simulations for the 
Heisenberg spin glass are unprecedented in system size. 
To achieve thermalization, we have needed not only a 
huge amount of CPU (7x 10 6 hours) but a careful samplc- 
by-sample thermalization check that allowed us to con- 
centrate efforts on the "hard" samples. The results for 



the spin-glass sector can be accounted for using only 
leading-order scaling corrections, but sublcading correc- 
tions are sizable for the chiral glass sector. This is the 
reason for the overestimate of Tqg m previous work i 10 ^ 1 
Data for L > 16 support the most economic scenario, 
Tsq = Tqg- We also see that the spin Binder parameter 
is not trivial at Tqg- 
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